% In this practical we will work with a single subject's EEG dat and perform an ERF
% analysis in sensor space. 
%

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

clear all; close all;

global OSLDIR;


%tilde='/Users/dantemant/Documents/meeg';
%osldir=[tilde '/osl1.2.beta.16dante'];    

tilde='/Users/woolrich/homedir';
osldir=[tilde '/matlab/osl1.2.beta.17'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

%testdir=[tilde '/meeg_data/meeg_signals/case_1084/111202'];
%testdir=['/Users/dantemant/Documents/meeg/duncan_data'];
testdir=[tilde '/vols_data/duncan_data'];

workingdir=[testdir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

clear spm_files_continuous spm_files_epoched;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[workingdir '/efMMspm8_EEG2_raw_0001.mat'];
spm_files_epoched{1}=[workingdir '/efMMspm8_EEG2_raw_0001.mat'];



% set up a list of mris
structural_files{1}='';

% defining experimental conditions and contrasts to be calculated

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to croat_settingse the (num_regressors x num_trials) design matrix:

fid_label.nasion='spmnas';        
fid_label.lpa='spmlpa';
fid_label.rpa='spmrpa';


%%%%%%%%%%%%%%%%%%%
%% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS
% Before running the beamformer we need to compute the forward model for
% each subject based on their structural scan.
% Make sure you check the results look reasonable!

for i=1:length(spm_files_continuous),

    %D=spm_eeg_load(spm_files_epoched{i});
    %fidnew=D.fiducials;
    %fidnew.fid.label

    S2=[];

    S2.fid_label=fid_label;        


    S2.D = spm_files_continuous{i};    % requires .mat extension    
    %S2.D = spm_files_epoched{i};    % requires .mat extension    
    
    S2.mri=structural_files{i}; % set S2.mri=''; if there is no structural available        
    S2.useheadshape=1;

    %S2.forward_meg='Single Shell';
    %S2.forward_meg='MEG Local Spheres';

    D=osl_forward_model(S2);

    %D=osl_neuromag_grad_baseline_correction(S2.D,'vector_view');

end;

%spm_eeg_inv_checkmeshes(D);
spm_eeg_inv_checkdatareg(D);
%spm_eeg_inv_checkforward(D, 1);



for ww=1:length(spm_files_epoched)

D=spm_eeg_load(spm_files_epoched{ww});

cnd=conditions(D);

mat=zeros(2,length(cnd));
for zz=1:length(cnd)
    str=cnd{zz};
    mat(1,zz)=str2num(str(5));   % here we assume that the format is LoadX_YY
    mat(2,zz)=str2num(str(7:end));
end

cat_var=unique(mat(1,:));

X=zeros(length(cnd),length(cat_var));

for k=1:length(cat_var)
    vect=find(mat(1,:)==cat_var(k));
    par_var=mat(2,vect);
    par_var=par_var-min(par_var)+1;
    X(vect,k)=par_var;
end

save([spm_files_epoched{ww} '.txt'],'X','-ascii');
design_matrix_summary={};
design_matrix_summary{1}=[spm_files_epoched{ww} '.txt'];


contrast={};

nc=length(cat_var);

for cont=1:nc
contrast{cont}=zeros(nc,1);
contrast{cont}(cont)=1; % main effects of a single condition (categorical variable)
contrast_name{cont}=['C' num2str(cont)];
end
for i=1:nc-1
    for j=i+1:nc
        cont=cont+1;
        contrast{cont}=zeros(nc,1);  % direct contrast between pairs of conditions
        contrast{cont}(i)=1; 
        contrast{cont}(j)=-1; 
        contrast_name{cont}=['C' num2str(i) '_vs_C' num2str(j)];
    end
end
for i=1:nc-1
    for j=i+1:nc
        cont=cont+1;
        contrast{cont}=zeros(nc,1);  % direct contrast between pairs of conditions
        contrast{cont}(i)=-1; 
        contrast{cont}(j)=1; 
        contrast_name{cont}=['C' num2str(j) '_vs_C' num2str(i)];
    end
end
cont=cont+1;
contrast{cont}=zeros(nc,1);  % average of all conditions
contrast{cont}(1:nc)=1; 
contrast_name{cont}='C_all';
% Currently there is only 1 subject.


% set up a list of SPM MEEG object file names (we only have one here)



% setup oat
oat=[];
%oat.source_recon.D_continuous=spm_files_continuous{ww};
oat.source_recon.conditions=unique(cnd);
oat.source_recon.D_epoched{1}=spm_files_epoched{ww}; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[]; % frequency range in Hz
oat.source_recon.time_range=[-0.1 0.5];
oat.source_recon.modalities={'EEG'};
oat.source_recon.method='beamform';
oat.source_recon.gridstep=12; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.fid_label=fid_label;
oat.source_recon.dirname=[spm_files_continuous{ww} '_wideband'];
%oat.source_recon.forward_meg='MEG Local Spheres';
oat.source_recon.forward_meg='Single Shell';
oat.source_recon.work_in_pca_subspace=0;
oat.source_recon.pca_dim=50;
% oat.source_recon.hmm_num_states=-1;
% oat.source_recon.hmm_num_starts=1;


oat.first_level.design_matrix_summary=design_matrix_summary;

for i=1:length(contrast)
    oat.first_level.contrast{i}=contrast{i};
    oat.first_level.contrast_name{i}=contrast_name{i};
end
oat.first_level.cope_type='acope';

oat = osl_check_oat(oat);  % function changed by DM


%% RUN THE OAT:

oat.to_do=[1 1 0 0];
oat.first_level.name=['wholebrain_first_level'];
oat.do_plots=1;
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% view the GLM design matrix
figure;imagesc(stats.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

% results = osl_get_recon_timecourses( osl_load_oat_results(oat,oat.source_recon.results_fnames{1});, mask_fname )

%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.time.reduce_time=1;
S2.resample_method='fft_ds';
S2.time_ds=2;
S2.time_range=[-0.1 0.5];
S2.data=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
[stats] = osl_reduce_data_to_visualize(S2);

S2=[];
S2.oat=oat;
S2.stats=stats;
S2.resamp_gridstep=12;
S2.first_level_contrasts=[1]; % list of first level contrasts to output
[statsdir,times,count]=osl_save_nii_stats(S2);

%% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.


% INSPECT THE RESULTS OF A CONTRAST IN FSLVIEW. Recall:
%S2.contrast{1}=[3 0 0 0]'; % motorbikes
%S2.contrast{2}=[0 1 1 1]'; % faces
%S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

con=1;
runcmd([getenv('FSLDIR') '/bin/fslview ' [statsdir '/tstat' num2str(con) '_12mm']  ' &']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end